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The general idea of a stochastic gauge representation is intro- 
duced and compared with more traditional phase-space expansions, 
like the Wigner expansion. Stochastic gauges can be used to obtain 
an infinite class of positive-definite stochastic time-evolution equa- 
tions, equivalent to master equations, for many systems including 
quantum time-evolution. The method is illustrated with a variety 
of simple examples ranging from astrophysical molecular hydrogen 
production, through to the topical problem of Bose-Einstein conden- 
sation in an optical trap and the resulting quantum dynamics. 

1 Introduction 

The original goal of pioneering physicists like Galileo, Newton and Einstein was 
to predict dynamical behavior in the universe - from projectiles to planets and 
even photons. However, complex systems in nature are not soluble with analytic 
or direct computational methods if the space of system descriptions is too large. 
Examples are master equations - widely used in many disciplines - or quantum 
theory itself, where the Hilbert space for many-body systems is enormous. This 
is one of the central problems of modern theoretical physics. 

The problem of complexity inspired a claim by Feynman^] that reads; 

• "Can a quantum system be probabilistically simulated by a classical uni- 
versal computer? ...the answer is certainly. No! " (Richard P. Feynman, 
Simulating Physics with Computers, 1982) 

In this paper, we will give a very general overview of new techniques for proba- 
bilistic simulation, called stochastic gauge representations[2j, which allow these 
'impossible' simulations. Like Wigner 's original representationjH], these are 
phase-space representations, but in phase-spaces of larger than classical dimen- 
sion. The purpose of the paper is to explain the abstract ideas behind stochastic 
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gauges, which are essentially an equivalence class of curved-space path-integrals 
in a complex phase-space. We compare the method with other phase-space 
approaches, as well as giving elementary examples. 

Even in Feynman's day, the introduction of Quantum Monte Carlo (QMC) 0] 
and related path integral techniques indicated a possible way to solve problems 
of large Hilbert space dimension. While these methods are used to treat canon- 
ical ensembles, they cannot treat real-time dynamics. Presumably, it was the 
failure of QMC methods in real time - due to the highly oscillatory behaviour 
of the real-time Feynman path-integral - that caused Feynman to neglect these 
types of random sampling methods. 

In a more recent publication, the time-domain quantum simulation problem 
was restated by Ceperley[n|, as follows: 

• "There are serious problems in calculating the dynamics of quantum sys- 
tems" (David M. Ceperley, Microscopic simulations in physics, 1999). 

The problems are due to the astronomical size of many-body Hilbert spaces. 
This makes it difficult to treat the quantum dynamics of Bose-Einstein con- 
densates - which typically have 10000 or more particles, and 10^'^'''^'' states in 
Hilbert space. In this paper, we show that stochastic gauge methods are versa- 
tile enough to treat complex quantum systems including master equations and 
canonical ensembles, as well as many-body quantum dynamics in real time. 

2 Phase-space representations 

Phase-space mappings, which map the discrete states of quantum theory into a 
classical-like phase-space, were originally introduced by Wigner in the form of 
the famous Wigner representation[3j. In the stochastic version of these methods, 
ensemble averages are mapped into trajectory averages - which can be numer- 
ically simulated. Phase-space representations have developed in three distinct 
stages. In the first stage, the Wigner representation, Husimij^ Q-representation 
and Glauber-SudarshanjJI |H] P-representation all use a classical phase-space of 
2M real dimensions for quantum systems corresponding to M classical modes. 
These methods typically will not give locally positive propagators for nonlin- 
ear quantum systems. Hence they cannot have a directly equivalent stochastic 
process. 

In the second stage, higher dimensional representations were developed. These 
include the Glauber R-representation[7| (which is non-positive), the Poisson rep- 
resentation (for incoherent master equations only), and for quantum systems, 
the positive P-representation[n|. The last two methods give locally positive 
propagators, by virtue of using a non-classical phase-space. However, while 
they work very well for damped systems, in highly nonlinear simulations they 
typically develop unstable trajectories ^1 leading to large sampling errors 
or even systematic errors in nonlinear quantum simulations. These simula- 
tion problems are caused by the distributions having power-law tails in phase- 
space^SI- When this happens, the distributions are not sufficiently localized to 
allow integration by parts, which is a crucial step in the derivation. 
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Phase-space 
Repn. 


Non- 
singular? 


Pos. 
dist? 


UV 
conv? 


2nd 
deriv? 


Pos. 
def? 


Stable? 


Wigner 


Yes 


No 


No 


No 


No 




Q 


Yes 


Yes 


No 


Yes 


No 




p 


No 


No 


Yes 


Yes 


No 




R 


Yes 


No 


Yes 


Yes 


No 




Pos. P 


Yes 


Yes 


Yes 


Yes 


Yes 


Indefinite 


Gauge 


Yes 


Yes 


Yes 


Yes 


Yes 


Yes 



Table 1: Table of representation properties required for stochastic quantum dy- 
namics 

In the current stage, gauge representations are used, which add a further 
complex gauge amplitude 0, to the representation phase-space. These provide 
a way to overcome Feynman's complexity dilemma, and are the subject of the 
present paper. For appropriate gauge choices these methods are conjectured 
to be exact in principle. They have no instabilities, and we will show that 
numerical simulations give no systematic errors, even in cases where previous 
methods would have boundary term errors. Numerical implementations are 
still necessarily approximate due to the usual limitations of finite computers, 
but error levels can be estimated and reduced to any desired level compatible 
with the hardware and time available. 

To understand the reason for these developments, it is useful to enumerate the 
requirements for a phase-space representation that can be used to map quan- 
tum dynamics into a stochastic differential equation, which can be simulated 
numerically. They are as follows: 

1. Non- singular: Essential for finite probabilities, although an initial delta 
function can be tolerated. 

2. Positive distribution: Needed to get positive initial probabilities. 

3. UV convergent: To control sampling error on lattices, vacuum fluctua- 
tions should not diverge at large momentum cut-off. 

4. 2nd-order derivatives: The mapping must generate at most second- 
order derivatives, to obtain diffusive phase-space behavior. 

5. Positive-definite propagator: The short-time propagator must be positive- 
definite, for an equivalent stochastic process to exist. 

6. Stable: Trajectories in phase-space should be stable to prevent boundary 
term errors: further restrictions on noise growth are also needed. 

How do the known phase-space representations compare? This is shown in 
Table 1, for a system of much current interest, the anharmonic oscillator — See 
SecO 
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It can be easily seen that the earlier phase-space techniques using a classical 
phase-space had other problems as well as not generating positive propagators. 
While the positive-P method removes almost all the problems of earlier tech- 
niques, it still has a disadvantage in that it can generate moving singularities 
from unstable phase-space trajectories. This is known to cause problems with 
large sampling errors, and boundarv-termspi H ITTl [T^ causing systematic errors 
in simulations of systems with extremely low damping. 

The gauge representation method[2], which unifies and extends some other 
closely related approaches fTK j [TH j I14|. removes this last problem by stabiliz- 
ing phase-space trajectories. It also points the way to the development of other 
phase-space representations in future, since the fundamental idea relies on very 
general properties of completeness, analyticity, and the existence of operator 
mappings into a second-order differential equation. 

If we extend our table to other cases, we see that some phase space rep- 
resentations appear more suited to calculations other than stochastic simula- 
tions. For example the symplectic tomography scheme of Mancini, Man'ko, and 
Tombesi^I], which expresses the quantum state as a probability distribution 
of a quadrature observable depending on a range of lab parameters, has been 
used to investigate quantum entanglement and failure of local realism, but has 
not to our knowledge led to many-mode quantum simulations, presumably due 
to the lack of a positive propagator in nonlinear evolution. The complex P 
representation[3 allows one to derive exact results for certain problems, but 
does not lead to stochastic equations, since the distribution is neither real nor 
positive. 

3 Stochastic gauges 

The idea of stochastic gauges can be summarized for the generic case of any 
linear time-evolution problem which, like quantum mechanics, can be expressed 
using a matrix product over a space of denumerable dimension. We first introduce 
a continuous basis Ao(a;) in the underlying Hilbert space of operators, with unit 
trace: Tr[Ao(Q;)] = 1. This must be an analytic function of the complex phase- 
space variables a, and have a mapping from the time-evolution problem that 
generates only second-order phase-space derivatives. The basic mathematical 
steps are as follows: 

1. Wish to solve dp/dt = 1^ p{t) where L is a matrix and p is a vec- 
tor of probabilities or amplitudes for the d distinct occupation numbers 
(A^i, . . . ,Nd), so its elements may be labeled pNi,...,Na- 

2. Introduce a (d -|- 1) dimensional complex phase-space: q = (0,a), with 
a renormalised (gauge) basis of analytic vector functions A(q) = r2Ao(Q:) 

3. Expand: p= J G{a)A{a)£('^+'^^ a 

4. Equivalent time-evolution using second-order derivatives: 
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5. Add arbitrary diffusion gauge vectors f (a) and drift gauges g(a) to give 
stability. 

6. Equivalent stochastic equation: da/dt = A. + 'B:<^{t) 
3.1 Ladder Operators 

To explain the procedure in more detail, consider a generic equation, which is 
typically a type of master equation for a quantum density matrix defined over 
a basis of number state occupation numbers: 

|p(0 = Lp(i) (1) 

where p has elements PNi,...,Nd, with the labels Nj corresponding to an occupa- 
tion of Nj in mode j. If the density matrix has off-diagonal entries, these can be 
regarded as elements of an enlarged vector, with d = 2M occupation numbers 
required for each entry in the case of M modes. Note that this linear problem 
is soluble in principle using diagonalization of L, but the (typically) large size 
of the matrix makes this impractical. 

Suppose L can be constructed from sums of products of matrix raising and 

lowering ('ladder') operators. These either increase (Lt) or decrease (L~) the 

-t- 

number of particles and multiply the probability by a function fj {Nj). Thus, 
for master equations, one might have: 



^^"/fe^.v-. * P) 



This general structure occurs in many physical systems, including Pauli-type 
master equations for positive probabilities found in many situations ranging 
from genetics to kinetic theory. In quantum problems, p is a density matrix in 
a number-state representation, while the ladder operators are bosonic creation 
and annihilation operators. The density matrix can be written as an enlarged 
vector on a basis of number-state projectors with d = 2M for the case of bosonic 
many-body theory on M classical modes. In this case, p(t) is not positive valued, 
but the time-evolution problem is still linear. 

There are also equations of identical structure found in 'imaginary time', 
which allow the calculation of canonical ensembles in many-body theory. In 
these cases, the operator norm is not preserved, but the equations are still linear. 
To give a simple example of the type of basis set that is of interest, a complete, 
analytic coherent state basis on the phase-space a = ( , a) = ( O , a, /?) for 
the case of a single harmonic oscillator is given by: 

A(a) = ^^|a)(r|e-"^ (3) 
where |q;) is a bosonic coherent state. 

3.2 Identities 

Identities can now be constructed that depend on the nature of the continuous 
basis set Ao{a). We only require that this basis set is an analytic function of the 
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continuous variables a. While it is common to use either Glauber coherent state 
projectors or Poisson distributions for this purpose, this is certainly not essential. 
More general basis sets like SU (N) coherent state projectors or general Gaussian 
bases are very likely to give even better results, as they often more closely 
approximate the physical quantum states of interest. 

Clearly, both raising and lowering identities are usually needed - in the follow- 
ing list, we indicate how the identities map matrices onto differential operators in 
the phase space. This is just a generalised version of the well-known equivalence 
between Heisenberg's matrix mechanics and Schroedinger's wave equation: 



LTA(a) 
L+A(a) 
A(a) 



Cjid,a) [A(«)] 
£+(a,a) [A(a)] 
ndn [A(a)] 



(4) 



Provided Ao(a) is analytic in a, one can use d = (Oq, d) to symbolize either 



d/dxj 



or 



d/dyj 



where a, 



+ iyj for j = 0, 1, . . . , d, and 



Xj as well as yj are real. These identities will be used later to specify which 
form of the derivative will be used to obtain a positive-definite diffusion term. 



3.3 Diffusion gauge 

Using the identities to eliminate ladder operators, we obtain an evolution equa- 
tion in integro-differential form: 

|p(t) = j G{cx)Ca [A(a)] d2('^+i)a . (5) 

Here the differential operator acts on the basis set, and must be of no more 
than second order: 

Ca = u + j^jd, + ]^Di,didj , (6) 

where the implicit summation is over i,j = l,...,d. At this stage, we notice 
that if we integrate by parts, we would obtain a possible solution to the time- 
evolution provided that boundary terms vanish, and that: 

^G(«,t) = £^[G(«)] . (7) 

Here the normally ordered differential operator is defined as: 

CN = U-djA'j + ^d.,djD,, . (8) 

This type of generalised Fokker-Planck equation is known to be equivalent to 
a curved-space path integral^Hl on the complex phase-space. This means that 
we have indeed reduced the dimensionality of the problem, in the sense that 
the phase-space dimensionality is much smaller than the dimensionality of the 
original vector. However, path-integrals are not always convenient for numerical 
calculations, and we wish to transform the equations further into a stochastic 
form. 
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As the first step toward a stochastic reformulation, define a d x d' complex 
matrix square root B called the noise matrix, where: 



D = BB^ 



(9) 



Since this is non-unique, one can introduce diffusion gauges^! [2] from a set 
of matrix transformations U[f(a)] with UU"^ = I. Using these, it is clear that 
given an initial square root B', it is always possible to construct another square 
root B, corresponding to an alternative 'diffusion gauge', with: 



B = B'U[f] 



(10) 



This is also an equally valid matrix square root, but it will in general have 
different stochastic properties. 

The matrices U[f] are just the usual set of complex orthogonal matrices where 
U[f]U^[f] = 1. These can either be fixed or variable functions of the phase- 
space coordinates, provided they satisfy growth restrictions. If we assume that 
they are d' x d' square matrices, then they are generated by the antisymmetric 
d' X d' matrices, that is, there are d' x {d' — l)/2 diffusion gauges. Generally, 
d' < d, since the diffusion matrix can have zero eigenvalues, although one can 
define larger noise matrices. 



3.4 Drift gauge 

A drift gauge Fokker-Planck equation is obtained as follows. Introduce d' arbi- 
trary complex functions g = {gj{a.,t) ), to give a new differential operator: 



(11) 



Here, Cq is equivalent to Ca, from the last identity in Eq Q. Summing 
indices over i,j = 0, . . . ,(i (where i,j = label the variable $7) , this can be 
rewritten in the form: 



G 



+ \Dijdidj 



(12) 



This removes the non-stochastic term U , and - with the correct choice of gauge 
- stabilizes the drift equations. 



4 Stochastic gauge equations 



Since the initial drift vector was A', the total complex drift vector, including 
gauge corrections is ^ = (C/il,A), where: 



A'-Bg 



(13) 



The total complex diffusion matrix ^ is a {d+1) x {d+l) matrix, with a new 
(d -t- 1) X d' square root B: 



D 





17g^B^ " 






BgJ7, 


BB^ 




B 



[fig,B^] =Rg 



(14) 
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4.1 Dimension-doubling 

We now introduce the technique used to produce positive-definite diffusion, 
which depends on the analyticity of the basis and the associated differential iden- 
tities. This technique is identical to that used for the positive-P representationf^, 
but is now extended to include the gauge amplitude variable as well. Define a 
2{d + 1) dimensional real phase space (xq, ..Xd^yo.. ,yd), with derivatives d^, 
where fj, labels all the 2{d + 1) real variables x^. 

Let: ^ = + and A = + iA^ with all the ^ and y forms real. 
Choose the alternative analytic forms of the differential operator so that: 

A^dj^£jd^+£.d] , (15) 



(16) 



With this identification of real derivatives, the original gauge differential oper- 
ator is written: 



(17) 



Next, on partial integration of the integral equation of motion, at least one 
valid solution for G must satisfy: 



d_ 
di 



G 



G 



By construction, the real diffusion matrix is a square of form: 



V 



BB^ 



(18) 

(19) 
(20) 



Clearly, P is positive definite. Hence, from the theory of stochastic equations [20] , 
provided some restrictions on growth are satisfied, one obtains the Ito stochastic 
differential equations: 

j^x^ = A^, + B^jCj{t), (21) 

where the real, Gaussian noise terms Cj(0 (for j = l,...,(i') are delta- 
correlated: 



m)Ut')) = Kt-t')5,j . (22) 
4.2 Central result of stochastic gauge theory 

Another, clearer way to write this result is to return to a complex vector nota- 
tion. Ito stochastic equations for the complex trajectory and gauge amplitude 
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17 are therefore obtained as follows: 



da 
IE 



17 [[/ + g • C(t)] 



A' + B:[C(t)-g] 



(23) 



Here, the arbitrary gauge terms g can be used to eliminate moving singu- 
larities that might be already present with the analytic drift A'; the actual 
simulated drift is A = A' — B : g. Gauges can be chosen freely to optimize 
simulations, in either real or imaginary time. 

It is essential to recognize that only the basis set, not the gauges, must be 
analytic functions on the phase space. Thus, while the original drift is usu- 
ally analytic, the gauge modified drift is best chosen not to be analytic. This 
is because the analytic continuation of systems of (generically) non-integrable 
nonlinear equations typically has moving singularities, which are trajectories 
that can deterministically reach infinity in a finite time. This is related to the 
Painleve conjecture of mathematical physics pH]- To remove the singularities, 
with their resulting boundary terms, a non-analytic gauge is therefore needed. 
Provided there are no boundary terms, all gauges are the same physically - but 
in practice, they give rise to different sampling errors, and therefore must be 
optimized carefully. 

Although removal of moving singularities appears necessary, this is not suffi- 
cient either to minimize sampling error or to guarantee the absence of boundary 
terms. In general, gauges must still be checked on a case by case basis. We have 
found that diffusion gauges that have no radial growth in the extended phase- 
space, and drift gauges where all phase-space trajectories are directed toward 
the origin at large enough radius, appear to eliminate boundary terms in most 
physically sensible examples. Clearly, a more rigorous investigation into these 
issues is still needed, since there may be anomalies even with these restrictions. 

4.3 Observables 

We typically wish to calculate physically observable quantities or moments of 
p in the form of a trace, where for this purpose the quantum density operator 
should be regarded as a matrix: 



If the problem involves Bose operators, with a coherent state basis, then 
the Hermitian observables can be written in a normally ordered form, O = 
OAr(a, a^^). The equivalent c-number expression for ensemble averaging is: 



(O) 



Tr[Op] 
Tr[p] 

/G(a)Tr[OA(Q)]d2(rf+i) 



cx 



0{a) = nON{a,P) 



(24) 
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so the quantum ensemble average ^Oyhas an immediate expression as: 

(a)=/^0(a)<i"'"s = ». (25) 

Here N = (O) g , to preserve the trace of the normahzed density matrix, and 
{O) g represents a stochastic average on the phase-space of all trajectories q , 
including the weighting factor f2 at each point in the trajectory. 



5 Examples 

The focus of this section is to give examples in simple cases that are exactly 
soluble, yet with statistics that are far from Poissonian. The reason for this is to 
illustrate how the stochastic gauge technique is successful in treating cases that 
would not be soluble using any previous simulation method, since the statistics 
are quite different to those of the underlying basis set. As examples, we will 
start with a simple chemical reaction master equation in which there are no 
quantum coherences, then move to a canonical ensemble example, and finally a 
quantum dynamics problem. 

None of the examples presents any real difficulties, since they are exactly 
soluble. However, our purpose here is to show that the stochastic gauge method 
gives correct results in cases where the solutions are already known. This, of 
course, is an essential first step toward treating more complex systems where the 
results are not known a priori. It is also interesting to see how these techniques 
have rather general applicability in physics and related scientific fields. This 
allows the possibility, for example, of combining imaginary time propagation 
for the initialization of the quantum system in a thermal ensemble, followed by 
real time propagation to simulate the response to a change in the Hamiltonian. 
Such experiments are common in many-body physics. 



5.1 Master equations 

The first example will be a type of Pauli master equation, in which there are 
only diagonal terms in the density matrix - so there are only simple proba- 
bilities in the original equations. These types of equation commonly occur in 
chemicalj22 and biological^22jj applications. A typical example is the astrophys- 
ically important problem of the formation of molecular hydrogen on interstellar 
grain surfaces jSBl- A simplified reaction model is then: 

2H -^k H2 
H -^^ H* 

This describes adsorption of hydrogen atoms via a rate (r) from an input flux 
H^^^^ , together with desorption at a rate 7. In addition molecule formation 
occurs at a rate of k. The corresponding master equation can be transformed to 
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Fokker-Planck form using the Poisson representation, giving an analytic steady- 
state solution for the m— th moment of n (the number of H atoms), in terms of 
Bessel functions: 

/■{0+) 
J — oo 

= (^) /^/fc+m-i(4xA72^)//vfc-i(4072fc) . (26) 
This gives the H2 production rate via = k (n^) . 
5.1.1 Stochastic equations 

Here we use an exact expansion of the distribution vector p using 'prototype' 
solutions, namely the complex Poisson distribution Ao(a): 

d 

[M<^)]m,...,N, = n K)^^- /N,\ (27) 
j=i 

By comparison, the original Poisson representation[2^ of Gardiner expands the 
distribution vector with a positive distribution of Poissonians, /(a), defined 
over a complex d-dimensional phase-space of variables ct, just as we do here 
- but without the extra weight-factor O. Including the weight factor, we can 
take advantage of the more general stochastic gauge procedure summarized 
in Eq Together with the corresponding Poisson identities, one finds the 

following Ito stochastic equations for molecule production including the gauge 
terms, where {C{t)C{t')) = ^{t ~ ''"')■ 



ngc 

[r - 7n - 2kn^] + iV2kn[C - g] . (28) 

If there is no gauge, the result of Fig (1) is obtained, corresponding to the 
original Poisson representation method. 

This result is clearly extremely inaccurate. It has a large sampling error in 
(n^) and we will see that it also has a systematic error in (n). The reason for this 
is that the original ungauged equations have an instability as n — > —00, leading 
to a moving singularity. This causes power-law tails and systematic boundary 
term errors in the resulting phase-space distribution, when there are small 
ratios. 

Fortunately, it is simple to stabilize these equations by adding non-analytic 
corrections to the drift. The simplest case is the 'circular' gauge, which replaces 
the analytic variable n by its modulus |n|: 

gc = iV2k{n — |n|) 



dil 

Ih 

dn 

d^ 
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Figure 1: Sampled moments of (n) (upper curve) and (n^) (lower curve) for 
astrophysical hydrogen molecule production in the Poisson represen- 
tation, parameters as in text. Adjacent hnes give upper and lower 
error bars caused by samphng error. The numerical values used here 
were k = 0.5, 7 = 0.1, r = 0.1 . 



Moment 


Analytic 


Poisson 


Circular gauge 


{Q) 


1.0 


1.0 


.994(10) 


in) 


0.407.. 


0.457(4) 


0.402(5) 


in') 


0.059.. 


0.094(8) 


0.061(2) 



Table 2: Table of observed moments, comparing analytic results with those for 
the Poisson representation and for the 'circular' stochastic gauge. The 
moment (n^) is critical for molecule production. Sampling error in 
brackets. 

In this gauge, the Ito equations are: 

^ = in{n - \n\) V2k C 

^ = r-n[j + 2k\n\]+inV2kC . (29) 

The corresponding results are shown in Fig (2), indicating a dramatically im- 
proved sampling error, and no systematic errors. 

For the circular gauge and for the Poisson expansion, the observed moment 
and its sampling error is given in Table 2, which tabulates the final near- 
equilibrium simulation results at t = 40, and compares them to the equilibrium 
analytic result for t = 00. For the stable circular gauge, the results are within 
one standard deviation of the analytic calculation in all cases. Other gauges 
are also possible - in fact, almost any gauge which suppresses the moving sin- 
gularities will give acceptable results, as long as no new boundary terms are 
introduced by the gauge itself. 

By comparison, the unstable ungauged Poisson method clearly gives enormous 
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<n>, <m> 
0.5 ^ 

0.4 - 

0.3 - 

0.2 - 

0.1 - 





8 16 24 32 40 

t 

Figure 2: Sampled moments of (n) (upper curve) and {n?) (lower curve) for 
astrophysical hydrogen molecule production in the 'circular' gauge, 
parameters as in text. Adjacent lines give upper and lower error bars 
caused by sampling error, which is invisible on this scale. 

sampling errors, with incorrect averages in (n), due to systematic boundary term 
errors. While this problem is relatively simple (for purposes of illustration), the 
stochastic techniques given here are easily extended to more complex problems 
where the original master equations cannot be solved directly. More details of 
this will be given elsewhere ^24j. 

5.2 Canonical ensemble 

For computational purposes, we can reduce the Bose gas Hamiltonian to a lattice 
Hamiltonian which contains all the essential features. This includes nonlinear 
interactions at each of M sites or modes, together with linear interactions that 
couple different sites together. Such problems are important for quantum gases 
trapped in optical lattices, or in low-dimensional environments, where evidence 
for quantum coherence and particle antibunching has been inferred in recent 
experiments. 

The simplest case that can represent a Bose-Einstein condensate (BEC) in a 
one-mode trap has M = 1, so: 

H = nw:n: +hx : : . (30) 

In this normally ordered Hamiltonian, the operator n = a^a is the boson 
number operator. The above Hamiltonian can be easily generalized to many 
important interacting Bose gas models. The canonical ensemble in thermal 
equilibrium for the one-mode case is an exactly soluble problem, which can be 
used to illustrate the gauge method. Applications in less trivial cases will appear 
elsewhere. It is an interesting historical note that the quantum correction to a 
classical canonical ensemble calculation was the first application jS] of the Wigner 
distribution, and Wigner regarded the effect of Bose or Fermi statistics to be a 
serious issue to be addressed in future. 
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The un-normalized grand canonical quantum density matrix is defined for our 
purposes as a slight modification to the usual form in statistical mechanics. We 
let: 



p = exp 



-tK - eN 



(31) 



where r = h/ksT, and K(^,a,a^) = H/h — p,N. We choose e <C 1 to give 
a high-temperature initial state at r = 0, with an initial occupation of no = 
l/[exp(e) — 1] « 1/e at each site. Thus, the effective chemical potential is 



K,N 



actually ^Ucfr = ^ — e/r. Since 
shows that the density matrix satisfies the equation 

d 



0, direct differentiation of Eq pT|l 



d- 



1 

2 L 



K,p 



(32) 



Solving this equation with the initial conditions at r = gives the solution for 
p{t) at lower temperatures, where quantum effects like Bose condensation will 



occur. 



Let us expand the density matrix p on an off-diagonal coherent state basis set 
in the manner of the positive P distribution. This is given in The initial 
G-distribution is Gaussian: 



Go (a) oc exp 



/no 6'^{a- P*)5^in-1) 



(33) 



To determine the effects of the 'Kamiltonian' K on G{a), it is first neces- 
sary to calculate the effect of the annihilation and creation operators on the 
projectors A(a). This is obtained as follows: 



a A (a) 
a^A(a) 
A (a) 



a A (a) 

[da + P] A{a) 



(34) 



together with the corresponding identities for the reversed orderings. Using 
these operator identities, the operator equation Ip^ can be transformed to a 
differential equation. The (ungauged) differential operator acting on the basis 
A (a) is 



[K{p, a, da + P)+ K{p, dp + a, (5)] 



-K{p,a,(3) + Y^ 



A'.dj + -D,d] 



(35) 
(36) 



To simplify notation, define n = a(3, ai = a, and 02 = P- Then A'- = (p 
2xn — ijj)aj/2 , and Dj 



Following the main procedure summarized in Eq (|23|1 . a gauge correction is 
utilized to stabilize coherent state paths in highly non-classical regions in phase- 
space. This allows one to benefit greatly from the over-completeness of coherent 
states, in reducing the sampling error and eliminating boundary terms. To 
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stabilize large modulus trajectories which otherwise can lead to boundary term 
errors and large sampling uncertainties, we choose drift gauges gj = i.^{n—\n\), 
giving: 



daj 
Ik 

dn 

d7 



-uj) /2-x \n\] aj + iaj^Cj 

n (37) 



2 



There is an intuitive physical interpretation. Since P = a* in the initial 
thermal state, each amplitude initially obeys a Gross-Pitaevskii equation in 
imaginary time, with quantum phase-noise due to the interactions. This causes 
non-classical statistics with a ^ (5* to emerge at low temperatures. Along each 
path an additional ensemble weight is accumulated, which is logarithmically 
proportional to the Kamiltonian K{a). The zero-temperature steady-state is 
the usual Gross-Pitaevskii approximation, together with quantum corrections. 

To illustrate the method, first consider the non-interacting case with x = 0) 
where we can set the gauge to zero, and define a = /3* , giving a diagonal 
Glauber P-distribution. Then: 

dn , , 

- = (^-^)n 

d^ ^, . , ^ 

— = n{fj,-uj)n (38) 
dr 

The power of the normal-ordered coherent state expansion is shown by the 
fact that these equations are deterministic, even though they include all quantum 
fluctuations. By contrast, the corresponding path-integral equations have large 
vacuum noise terms. These equations can be integrated immediately to give a 
Bose-Einstein thermal ensemble with Gaussian fluctuations : 

^ {n{T))s exp([u;-MefTH-l ■ ^ ' 

Next, consider the exactly soluble interacting case, involving a single Bose 
mode with: 

H{a,a^) =hx:n^ : ■ (40) 

A numerical simulation is of most interest here, as it can be generalized to 
other Bose gas systems of greater complexity. It is straightforward to obtain 
agreement with the exact solution for large boson number, as quantum-noise 
corrections are small in this limit. Instead, we focus on the case furthest from co- 
herent statistics with /J, = X = 0-5, giving just one boson in the zero-temperature 
limit, and choose e = 0.1. This case is shown in Fig H5.2p . as well as a comparison 
with the exact results. 



The results can be seen to agree well with the exact analytic ones, and the 
agreement is easily improved by increasing the number of stochastic trajectories. 
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Figure 3: Simulation (solid line) versus exact results (dotted line) for the boson 
number n = (n) and correlation function g2 = {-.v? :) /v? of the exactly 
soluble anharmonic oscillator case with parameters X = A* = 0-5, £ = 
0.1. 

This excellent agreement also occurs with much larger numbers of bosons. The 
physical behaviour of strong particle antibunching is in agreement with evidence 
deduced from recent experiments in optical lattices. 

5.3 Anharmonic oscillator 

The quantum dynamics of the anharmonic oscillator is the subject of much 
current interest. One can combine the previous canonical technique with a real- 
time evolution, in order to model a transient experiment in which a BEG is first 
cooled, then allowed to evolve after a change in the Hamiltonian. We consider 
the Hamiltonian to be Eq JHOI), the same as previously. 

This type of problem has been studied previously as part of more extended 
multi-mode studies on quantum solitons, leading to the prediction and obser- 
vation of quantum squeezing in solitons I26| . Other applications include 
first-principles studies of evaporative coolingj22], and a treatment of phase- 
diffusion|28| using an approximate Wigner technique, starting from a coherent 
state. The present technique is exact rather than approximate, though it pre- 
dicts, as expected, the same behaviour of phase-diffusion and amplitude decay 
from an initial coherent state. However, we shall just consider the one-mode 
case here. 

We find the following salient points: 

• Diffusion gauges work better than drift gauges at controlling sampling 
error 
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• Sampling error grows with time 



Diffusion gauges work by trading off (increasing) phase noise against in- 
tensity noise 

Even better results occur if drift and diffusion gauges are combined 
It is possible to simulate past the time of amplitude decay 



As previously, we use the stochastic gauge procedure of Eq l(23|l . The resulting 
Ito equations in an arbitrary drift gauge are as follows, where a, [5 are the two 
variables that correspond to a, a) in a coherent state (positive-P) expansion. We 
define r = xt-, w = 0, and n = a(3, corresponding to n = a^a: 



da 

Ih 

dp 

dn 



—2ina + (1 — i)a(Ci ~ 9i] 

2in(3 + (1 + i)/3(C2 - 92) 
2 



(41) 



Here (O(r)O(rO) = <5ii<5(r - r'). 

These equations correspond to a diagonal noise matrix of form: 



B' 



(l-i)Q, 
0, {l+i)(5 



(42) 



It is convenient to use an equivalent diffusion gauge with a noise matrix B = 
B'U defined in terms of a parameter A, using an orthogonal transform U so 
that: 



U 



cosh. A, —ismhA 
i sinh A, cosh A 



(43) 



We also introduce new variables 6, 6 where: 



n 



a/ 13 



o«0 



(44) 



These variables are interpreted as the logarithmic amplitude and phase re- 
spectively. Their equations including gauge terms are: 



d^ 
dB_ 

d^ 

dn 

d^ 



2 -in- e^(l + i) [(Ci - 91) - i{C2 - 92)] 



e-^{l-i)[{Ci-9i)+i{C2-92)\ 
2 



(45) 
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Figure 4: Graph of mean sampling error in expectation values of quadrature 
moments versus time, using several different types of gauge. The 
combined diffusion and drift gauges give the best result, with an order 
of magnitude longer time duration before the sampling error becomes 
substantially large. Coherent state initial conditions with amplitude 
a(0) = 3. Constant diffusion gauge A = 1.4, 10^ trajectories. 



Without any drift gauge, these equations have the problem that the loga- 
rithmic ratio of amplitudes (imaginary phase) can grow rapidly whenever 
the quantum noise causes n to have an imaginary part, resulting in a large 
sampling error. This can be controlled to some extent with the diffusion gauge 
technique [inj, which allows us to make A large and positive, thereby reducing 
the quantum noise in n. This, however, is at the expense of an increase in 
quantum noise in the phase, and can only delay the onset of the rapid growth 
in 

Even better results are obtained by combining diffusion and drift gauges. We 
choose the gauge gi = {1 + i)'is[n]e~^ = —ig2- This has the property that only 
noise drives the ratio of amplitudes 



dr 

as opposed to the un-gauged equation of: 



(C2 - Ci) (46) 



d9 



49[n]+e^(C2-Ci) (47) 



dr 



Typical results are shown in Figure Q, which gives the time-dependence of 
the sampling error with different gauges. This indicates substantial increases in 
useful simulation times compared to any previous phase-space technique |14[ll6| . 
We note that the sampling error still increases with time. This appears to be 
inevitable with current stochastic gauge methods, which stabilize trajectories 
but introduce an increasing uncertainty in the relevant quantum amplitude. 
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6 Summary 



Complexity is a major, even central problem in modern theoretical physics. Our 
fundamental description of nature, quantum field theory, is incredibly complex 
in terms of the Hilbert space dimension - far more so than any classical descrip- 
tion. Mappings to phase-spaces of reduced dimensionality therefore provide an 
attractive route to overcoming this complexity problem. For best results, how- 
ever, one wishes to have a stochastic mapping, where the dynamics in phase- 
space obey a local stochastic description. In these cases the mapping permits a 
way to sample the complex dynamics over a finite set of samples, thus providing 
a controllable approximation to the exact dynamics. 

The Wigner representation is the pioneering method in phase-space that maps 
quantum mechanics into a classical phase-space. However, it has some draw- 
backs. It gives large vacuum fluctuations in quantum field simulations, and is 
essentially non-stochastic, as it is not a positive-definite representation. All the 
other techniques that have been introduced for classical phase-spaces, like the 
Glauber P-representation, have similar problems. This includes even the Q- 
representation - which is always positive, but has no local positive propagator 
when the Hamiltonian is nonlinear. More modern techniques like the positive-P 
representation, solve most of the technical problems due to lack of positivity, but 
can give boundary term errors due to moving singularities in the drift equations. 

Gauge techniques solve these known mathematical problems by stabilizing the 
drift equations. However, in specific cases, one still needs to find the optimum 
method or gauge field. Several examples of workable gauges in different cases 
and bases have been given here, and their utility demonstrated by comparison 
with exact results. In the long run, such techniques can allow one to make 
progress toward calculations that involve interacting many-body systems. We 
note that the gauge approach is rather general, since in any particular case one 
can optimize the basis set, the gauge, and even the algorithm for sampling the 
weighted trajectories. 

Prom the larger theoretical perspective these are important issues, as we treat 
more challenging complex systems, including possible tests of quantum mechan- 
ics in new regions of macroscopic and entangled quantum systems. 
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